########################################################
#Figure2
#This figure plots the empirical cumulative distribution, separately for households in treatment and control blocks, of value received divided by value entitled per month, pooling the endline one months

#########################################################

#set colors
cColor <- "#97B8DE"
tColor <- "#4669BD"


#Read in EL data
ELdat <- readstata13::read.dta13(paste0(SurveyDataDir, "JH_ePOS_HH_DataforAnalysis.dta"), convert.factors=FALSE)

ELdat[ELdat == "."] = NA


############################
# Value received per month
############################
tmp.df_value <- as.data.frame(
  na.omit(as.data.frame (cbind(ELdat$value_total_jan17, 
                               ELdat$value_total_feb17,
                               ELdat$value_total_mar17,
                               ELdat$treatment,
                               ELdat$uid,
                               ELdat$ghost_final))))



names(tmp.df_value) <- c("value_jan17","value_feb17","value_mar17","treatment", "uid","ghost_final")



#Reshape to long
tmp.df_value <- reshape(tmp.df_value, varying = c("value_jan17","value_feb17","value_mar17"),
                       v.names = "value",
                       timevar = "month",
                       times = c("jan","feb","mar"),
                       idvar = c("uid","treatment"),  direction = "long")
tmp.df_value$value <- as.numeric(tmp.df_value$value)

################################
#Statutory entitlement value
################################

tmp.df_statentvalue <- as.data.frame(
  na.omit(as.data.frame (cbind(ELdat$stat_ent_value_total_jan17,  
                               ELdat$stat_ent_value_total_feb17,
                               ELdat$stat_ent_value_total_mar17,
                               ELdat$treatment,
                               ELdat$uid,
                               ELdat$ghost_final))))

names(tmp.df_statentvalue) <- c("stat_ent_value_jan17","stat_ent_value_feb17","stat_ent_value_mar17","treatment", "uid","ghost_final")


tmp.df_statentvalue <- reshape(tmp.df_statentvalue, varying = c("stat_ent_value_jan17","stat_ent_value_feb17","stat_ent_value_mar17"),
                              v.names = "stat_ent_value",
                              timevar = "month",
                              times = c("jan","feb","mar"),
                              idvar = c("uid","treatment"),  direction = "long")

tmp.df_statentvalue$stat_ent_value <- as.numeric(tmp.df_statentvalue$stat_ent_value)


tmp.final <- merge(tmp.df_value, tmp.df_statentvalue, by=c("uid","treatment","ghost_final","month"), all=TRUE)

#####gen proportion var ###
tmp.final$prop_stat <- tmp.final$value/tmp.final$stat_ent_value

####### store mean ######################
mean_stat <- c(mean(tmp.final$prop_stat[tmp.final$treatment == "0"],na.rm=T), 
               mean(tmp.final$prop_stat[tmp.final$treatment == "1"],na.rm=T))
means_stat <- as.data.frame(mean_stat) 
means_stat$treatment <- c(0, 1)

tmp.final <- tmp.final[!is.na(tmp.final$prop_stat),]

#check percentiles

quantile(tmp.final$prop_stat, c(.50, .85, .90, .95)) 


########### output graph ##############

ggplot(tmp.final) + 
  stat_ecdf(aes(x = prop_stat, colour = factor(treatment)),size=1.2) +
  ylab("Proportion") + xlab("Value delivered/statutory value entitled per month")  + 
  theme(plot.title = element_text(lineheight=.8, size=rel(1.5))) +
  scale_colour_manual(values = c(cColor, tColor), 
                      breaks=c(0,1), 
                      labels = c("Control", "Treatment"))  +
  theme(legend.key = element_blank(), panel.background = element_blank(),
        panel.grid.major = element_blank(), panel.grid.minor = 
          element_blank(), panel.border = element_blank(),
        legend.title=element_blank(), legend.position="bottom",
        legend.key.size = unit(1, "cm"), legend.key.width = unit(2, "cm"),
        legend.text = element_text(size = 17, colour = "black"),
        axis.title=element_text(size=17), text = element_text(size=17),
        axis.ticks.length=unit(0.3,"cm")) +
  geom_vline(data=means_stat, aes(xintercept=mean_stat, colour=factor(treatment)), alpha = 0.8,
             linetype="dashed", size=1.2) + coord_cartesian(xlim=c(0, 1.04)) + 
  scale_x_continuous(limits=c(0,1.04),breaks = seq(0,1,.2),expand=c(0,0))+scale_y_continuous(expand = c(0, 0)) 

#cut off at 90% percentile
ggsave(file=paste(OutputDir, "Figure2.pdf", sep=""), width=12, height=8)






